#!/bin/env Rscript

library(HMMcopy)

argss <- commandArgs(trailingOnly = TRUE)

print(argss)


samplew <- argss[2]
sample_name <- argss[3]
gcw <- argss[4]
mapw <- argss[5]

# Correct readcount and estimate copy number
sample_uncorrected_reads <- wigsToRangedData(samplew, gcw, mapw)
# control_uncorrected_reads <- wigsToRangedData(controlw, gcw, mapw)

# The correctReadcount requires at least about 1000 bins to work properly
# process with chr separately?
sample_corrected_copy <- correctReadcount(sample_uncorrected_reads)
# control_corrected_copy <- correctReadcount(control_uncorrected_reads)

# Segmentation

# Normalizing sample by control
# sample_corrected_copy$copy <- sample_corrected_copy$copy - 1
# Segmenting
segmented_copy <- HMMsegment(sample_corrected_copy)
# retrieve converged parameters via EM
# param <- HMMsegment(sample_corrected_copy, getparam = TRUE)
# param$mu <- log(c(1, 1.4, 2, 2.7, 3, 4.5) / 2, 2)
# param$m <- param$mu
# perform segmentation via Viterbi
# segmented_copy <- HMMsegment(sample_corrected_copy, param)

# log2 states to normal states 
sample_corrected_copy$copy=2^(sample_corrected_copy$copy+1)
segmented_copy$segs$median=2^(segmented_copy$segs$median+1)


# Visualization

chr <- c(1:22, "X", "Y")

source(paste(Sys.getenv('tools_path'), 'script/cnv_p_hmmcopy.r', sep='/'))
cnvphmmcopy(sample_corrected_copy, segmented_copy, chr=chr, sample_name=sample_name)

q(save = "yes", status = 0, runLast = FALSE)
